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We present a fully quantum many-body treatment of dark solitons formed by ultracold bosonic 
atoms in one-dimensional optical lattices. Using time-evolving block decimation to simulate the 
single-band Bose-Hubbard Hamiltonian, we consider the quantum dynamics of density and phase 
engineered dark solitons as well as the quantum evolution of mean- field dark solitons injected into 
the quantum model. The former approach directly models how one may create quantum entangled 
dark solitons in experiment. While we have already presented results regarding the latter approach 
elsewhere [Phys. Rev. Lett. 103, 140403 (2009)], we expand upon those results in this work. In 
both cases, quantum fluctuations cause the dark soliton to fill in and may induce an inelasticity in 
soliton-soliton collisions. Comparisons are made to the BogoUubov theory which predicts depletion 
into an anomalous mode that fills in the soliton. Our many-body treatment allows us to go beyond 
the Bogoliubov approximation and calculate explicitly the dynamics of the system's natural orbit als. 
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I. INTRODUCTION 



Solitons are robust nonlinear waves that appear in 
many facets of nature. They are localized, persistent 
in time, and collide elastically. Solitons with both phase 
and amplitude are often described to lowest order in sys- 
tems with weak local nonlinearity by the famous nonlin- 
ear Schrodinger equation (NLS) [1]. For example, the 
NLS describes twist and curvature of a string in three 
dimensions [2]. However, corrections beyond the NLS, 
necessary in order to more closely model soliton proper- 
ties, are less universal. For instance, for solitons in optical 
fibers one must consider higher order dispersion. Bose- 
Einstein condensates (BECs) are also well-known to be 
modeled by the NLS, or Gross-Pitaevskii equation (GP). 
Dark solitons were theoretically predicted in BECs 0, Q 
and subsequently have been observed in a number of ex- 
periments [i, i, 0, i, i, [TO, [H . The NLS then represents 
taking the lowest order average of quantum operators 
from the quantum many-body description; therefore, it is 
called a mean-field theory. One naturally expects thermal 
or quantum corrections to mean-field theory. Regar ding 
thermal corrections, recent finite temperature studies [12j 
based on coupling Boltzmann equations to the NLS have 
indicated that thermal fluctuations cause dark solitons 
to decay. The question then arises, are quantum correc- 
tions, called quantum fluctuations, sufficiently large so as 
to change the basic properties of dark solitons, a funda- 
mental prediction of mean-field theory and a ubiquitous 
property of weakly nonlinear systems? 

The answer so far has been a resounding no. First, 
BEC soliton experiments have worked in regimes with 
very small quantum fluctuations. Such fluctuations can 
be calculated quantitatively with the Bogoliubov-de- 
Gennes equations (BDG) for weakly interacting BECs, 



and one finds much less than 1% for the system's ground 
state [l3|. Second, BDG methods have been used to de- 
scribe quantum fluctuations of solitons. The main result 
is that within BDG framework the dark soliton, which 
manifests as a density notch with a phase jump across 
it, appears to be filled in: the soliton's position is un- 
certain, so that in an ensemble average of measurements 
it appears blurred [13, [lH, [TO] • However, these studies 
predict that in any particular measurement the soliton is 
indeed localized, and the predictions of mean-field the- 
ory for the density, phase, and stability properties of the 
soliton hold. 

In this article, we provide new answers to the question 
of how quantum fluctuations affect dark solitons, expand- 
ing on the results which two have us have presented in 
Ref. [13] • We show that in a regime of stronger quantum 
fluctuations, the results based on BDG theory no longer 
hold. Solitons are not just delocalized in an ensemble 
measurement, but do in fact decay in the zero temper- 
ature quantum setting. They even collide inelastically. 
By use of a recently invented method for entangled quan- 
tum many-body dynamics in one dimension [l^l , we show 
explicitly how both mean-field theory and BDG theory 
break down. To obtain large quantum fluctuations we 
make use of an optical lattice [H, [l^ . BECs in optical 
lattices are known to have two distinct quantum phases, 
the Mott insulator and the superfluid [21"] . We emphasize 
that throughout our study we remain in the superfluid 
phase. We also treat higher order correlations, a key 
signal for dark soliton decay in a fully quantum theory, 
where the mean field does not force symmetry breaking 
on the system [H, [H, [^ • 

In general, the unprecedented tunability permitted by 
optical lattice systems gives both theorists and experi- 
mentalists alike an ideal setting for investigating quan- 
tum many-body phenomena in a controlled fashion. In 
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FIG. 1: (Color online) Density and phase engineering a stand- 
ing soliton. To engineer a standing soliton in a lattice, we (a) 
take a free lattice, (b) apply a blue-detuned Gaussian beam 
to dig a density notch, (c) turn off the Gaussian beam and 
raise half the lattice to imprint a phase, and (d) lower the 
lattice to its original state. The lattice heights are the same 
in each case, the number of atoms indicates the density, and 
red indicates atoms lower in phase by amount tt. 

these systems, experimentalists can adjust the height of 
the lattice by changing the intensity of the laser field 
used to create it. With this knob, one can vary the ra- 
tio of the interaction to kinetic energy in the governing 
lattice Hamiltonian, e.g. the Bose-Hubbard Hamiltonian 
(BHH). In fact, the first major experimental result for 
BECs in optical lattices [25[ used just such a carefully 
controlled increase in optical lattice height past a crit- 
ical value to demonstrate the Mott-superfluid quantum 
phase transition in the BHH [2l|, [2^ . More generally, at 
present, experimentalists also have precise control over 
the lattice site filling factor, lattice geometry and di- 
mensionality, species type, and atom-atom interaction 
strength, symmetry, and sign; the latter properties are 
due to use of a Feshbach resonance. The system is also 
well-isolated from the environment, greatly curbing the 
intervention of finite- temperature effects. All of these ap- 
pealing properties make it possible to investigate out-of- 
equilibrium quantum many-body dynamics governed by 
the corresponding Hamiltonian, an aspect of the problem 
on which we focus in this paper and one that received lit- 
tle attention before the advent of optical lattices. Specif- 
ically, we study how quantum many-body phenomena, 
e.g., quantum fluctuations and quantum entanglement, 
affect the behavior of dark solitons formed by ultracold 
bosons loaded into one-dimensional optical lattices. 

Regardless of the physical system in which they oc- 
cur, solitons are characterized by stable, non-dispersing 
propagation of either a minima (dark solitons) or a max- 
ima (bright solitons) of the physically relevant depen- 
dent variable, e.g., particle number density in BECs. It 
has been well-known for decades, since the seminal pa- 
pers by Tsuzuki on nonlinear waves in the CP equa- 
tion [23] and by Zakharov and Shabat on the inverse 
scattering transform [2^, [2^ , that solitons exist as solu- 
tions to the one-dimensional (ID) nonlinear Schrodinger 
equation (NLS). It is also well-established that in the 
mean-field limit Bose-Einstein condensates (BECs) are 
well-described by the NLS [l^. The success of the GP 
equation as a model to describe the ultracold atomic Bose 
gas, as well as the observation of solitons in systems such 



as photonic crystals and nonlinear waveguide arrays, has 
excited a broad new interest in solitons in the last decade. 
Both dark [sl, @] and bright [sO, [3l| solitons were ob- 
served rather early on in systems of BECs in harmonic 
trap geometries, and, due to ever improving experimental 
technique, there have very recently been exciting new de- 
velopments regarding the creation and evolution of dark 
solitons formed by BECs [7, 8, 9]. 

To treat solitons in the context of a lattice via mean- 
field theory, there are two approaches: (1) a full solu- 
tion of the NLS with an external lattice potential [32], in 
which case the lattice soliton solutions of the NLS have 
been mapped out in detail [sl,!!^, or (2) a finite-element 
discretization of the NLS obtained by projecting the con- 
densate wave function onto a basis of localized wave func- 
tions at each lattice site Q thus resulting in the discrete 
nonlinear Schrodinger equation (DNLS). Discrete lattice 
solitons have been studied using the DNLS both in purely 
mathematical [36|, [s^] and BEC-related [35|, [s^ contexts. 
On the experimental side, bright gap solitons have been 
observed in a system o f rep ulsive ^^Rb atoms trapped in 
a weak optical lattice [3^]. We emphasize that we will 
treat dark solitons with a size much greater than the lat- 
tice spacing, not gap solitons deep in the band gap which 
are localized on a few sites. 

In contrast to these mean- field approaches, we use the 
BHH, a fully quantum model. The DNLS emerges as the 
equation of motion when taking the mean-field limit of 
the BHH [40]. This is why we are able to make quanti- 
tative predictions regarding the accuracy of using mean- 
field theory to describe soliton dynamics. Related studies 
regarding the effect of quantum fluctuations on excited 
condensates in continuous geometries, e.g., dark solitons, 
using few-mode approximations and BDG-based meth- 
ods can be found in the works of Dziarmaga, Sacha, and 
Karkuszewski [3, [H, EB, Sll, |42[, as we have already 
alluded to. Using related methods. Law studied the dy- 
namical depletion of dark solitons created via phase im- 
printing [43]. Below, we make connections with these 
works in our study. An early study of quantum lattice 
solitons is presented in |4J] ; a quantum theory of solitons 
in optical fibers was developed by Lai and Haus [45]. 

Thus, in conjunction with Ref. [l^, we present the 
first full entangled dynamical quantum many-body treat- 
ment of dark solitons formed by ultracold atoms in opti- 
cal lattices. In SecHIl we introduce the pertinent models 
used and briefly summarize the techniques used for solv- 
ing them. In Sec. Illli we define all of our measures for 
system characterization. In Sec. IIVI we present results 
of quantum evolution of dark solitons obtained by den- 
sity and phase engineering ground states of the BHH, a 
method in complete analogy to how dark solitons could 
be created in optical lattices experimentally. In Sec. 
we consider quantum evolution in the BHH of mean-field 
solitons generated in the DNLS and add new measures 
and insight to the results in [13] • In Sec. IVIi we per- 
form a Bogoliubov analysis of the dark soliton solutions 
in Sec.[Vl Finally, in Sec. lVIIi we summarize our results. 
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II. MODELS AND METHODS 



The governing model we employ throughout this paper 
is the single-band BHH [2l|, [2y] in one spatial dimension: 
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H = -J^(St^^6fe + h.c.) 

^ M M 
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(1) 



where J is the nearest-neighbor hopping coefficient, U is 
the on-site interaction coefficient, Ck is an external po- 
tential aside from the lattice potential, and M is the to- 
tal number of lattice sites. The destruction (creation) 
operator hk {h\) destroys (creates) a boson in the low- 
est Bloch-band Wannier state localized at lattice site /c, 
while the number operator hk = i\^k counts the boson 
occupation at site k. We employ box (open) boundary 
conditions in all calculations, hence the restriction on the 
hopping term to contain only M — 1 terms. 

It is straightforward to derive Eq. ([T|) from the contin- 
uous second-quantized many-body Hamiltonian describ- 
ing a system of identical bosons with binary interactions 
after making the following assumptions: (1) the bosons 
interact with a short-range contact delta function poten- 
tial; (2) the tight-binding approximation is valid so that 
next-nearest-neighbor hopping and nearest-neighbor in- 
teractions can be neglected; and (3) all relevant energy 
scales, e.g., temperature, hopping, and interaction ener- 
gies, are smaller than the band spacing so that we can 
restrict the model to include only the lowest band. This 
must be dynamically true as well as statically. 

From this derivation, one can write down the coeffi- 
cients J, [/, and Ci in terms of the single-particle Wannier 
functions on which Eq. ([T]) is based [26| . Equation ([T]) is 
a truncation of the full continuous Hamiltonian that very 
accurately describes ultracold bosons trapped in optical 
lattices. We note that the ID model is obtainable experi- 
mentally from a full 3D lattice by increasing the strength 
of the lattice potential in the transverse directions, thus 
suppressing hopping in those directions. This results in 
an array of ID tubes with a lattice potential in the longi- 
tudinal direction. This is not a true ID quantum model 
though, in that we assume the transverse confinement 
length is still much larger than the range of the inter- 
atomic interaction potential and the actual atomic size. 
That is, even though excitations can occur only along 
the longitudinal direction, the underlying scattering pro- 
cess is still three-dimensional, i.e., we are in the quasi- ID 
regime. 

As alluded to in Sec. [H the DNLS can be obtained 
either as a finite-element discretization on the lattice of 
the continuous NLS [35[ or as a mean-field approximation 
of the BHH [io^. In the latter case, it is assumed that the 
many-body state takes the form of a product of Glauber 



coherent states at each site k: 

M CO / 

l*) = (g)kfe), where j^^) ^ e'l^'^l'/^ V i£^|n). 

(2) 

Then, after taking the expectation value of Heisenberg's 
equation of motion, ihdtbk = [bk^H], with H given by 
Eq. ([1]) but in normal ordered form, we arrive at the 
DNLS [46]: 

ihipk = - J (V^/c+1 + V^fe-i) + U\ilJk\^i^k + efeV^fe, (3) 

where Zk = {H) = i^k is the time-dependent coherent 
state amplitude at site k that is equivalent to the conden- 
sate order parameter in the limit of a product of coherent 
states. In Sec.|Vl we use a truncated form of the coherent 
states in Eq. ([2]) to build quantum analogs in the BHH of 
the dark soliton solutions of the DNLS, while in Sec. IIVI 
we model our quantum soliton engineering in the BHH 
after an analogous calculation in the DNLS. Finally, in 
Sec. (Vll we treat the stability of DNLS soliton solutions 
using the Bogoliubov method. 

For all simulations of the BHH in both real and imagi- 
nary time, we employ the time-evolving block decimation 
(TEBD) algorithm 13, [13] that allows efficient classical 
simulation of the dynamics of slightly-entangled quantum 
states in ID lattice Hamiltonians. The main convergence 
parameter, denoted Xi ^^is numerical method is the 
number of basis sets retained in a Schmidt decomposi- 
tion at each bipartite splitting of the lattice. For the 
specifics of our implementation of the algorithm, we re- 
fer the reader to Refs. [4^,[4^. The other approximations 
inherent in the algorithm are (1) a truncation of the lo- 
cal Hilbert space dimension at d states so that each site 
can contain at most d—1 bosons and (2) a second-order 
Suzuki- Trotter representation of the time-evolution oper- 
ator e~^^^^^^ over small time steps 5t. In all results pre- 
sented herein, we have explicitly checked for convergence 
in x-> and St. On the other hand, for all simulations 
of the DNLS in both real and imaginary tinie, we use a 
standard fourth-order Runge-Kutta routine [4g . 



III. MEASURES 

We now clearly define the measures that we use to 
probe the physical properties of the system during time 
evolution. The average particle number (n^), i.e., the 
expectation value of the number operator at site /c, 
is a local observable which corresponds directly to the 
particle density observed in experiment averaged over 
many runs. Another local number-related observable 
is the normalized particle number variance defined as 
r]k = {{Auk)'^) / (hk) = {{hi) - {hk)'^)/{hk). This mea- 
sure characterizes the deviation of on-site number statis- 
tics away from the classical Poissonian limit. Specifically, 
?7/c = 1 for atom-number Glauber coherent states which 
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obey Poissonian number statistics, rjk < I for a number- 
squeezed state, 77/e > 1 for a phase- squeezed state, and 
r]k = for a single Fock state of definite occupancy. The 
order parameter (6/^), defined as the expectation value 
of the boson destruction operator at site /c, corresponds 
directly to the DNLS solution t/jk in the limit of a direct 
product of atom-number Glauber coherent states. 

Within our d-dimensional local Fock space, we de- 
fine a Hermitian phase operator, according to the Pegg- 
Barnett [49] prescription as 



one another so that the state of a given particle is, in gen- 
eral, entangled with the remaining N — 1 particles |5Q| |. 
We refer to this type of entanglement as particle entan- 
glement. For each type of entanglement, we consider the 
average impurity Q and the von Neumann entropy S'vn- 
For the case of spatial entanglement, we define the aver- 
age local impurity \E^\ as 
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(4) 



where we have set a reference phase Oq to zero. We use 
(^/c), the expectation value of the on-site phase operator, 
as a quantum measure of average relative phase between 
lattice sites in characterizing quantum solitons. 

The lattice representation of the single-particle den- 
sity matrix is an M x M matrix with elements defined as 
(Pgp*)ij = (b^jbi), which can be used to calculate the sys- 
tem's single-particle natural orbitals and quantum deple- 
tion. The (j + 1)*^ most highly occupied natural orbital 
is computed as the (j + 1)*^ largest eigenvector of p^^. 



where pk is the reduced density matrix describing site 
k. The measure Qmodes quantifies multipartite entangle- 
ment in the sense that it averages the bipartite entangle- 
ment between each site and the remaining sites. We note 
that the average local impurity is equally well described 
as the average local mixedness, a particular case of the 
generalized entropy [52I. [53l . [5^ , or a local average quan- 
tum Tsallis entropy [55|. On the other hand, the local 
von Neumann entropy defined as 



5vN,/e 



-Tr[p,log,(pfe)] e [0,1], 



(6) 



measures the entanglement between site k and all other 
sites. The analogous quantities for particle entanglement 
are the single-particle impurity, 



IS 



The k component of this eigenvector, denoted 
the coefficient of the natural orbital associated with site 
k; the eigenvalues of p^^ are the average occupation num- 
bers of the natural orbitals and are denoted by Nj with 
^0 ^ ^1 ^ • • • ^ Nm-1' We assume eigenvectors nor- 
malized to unity: Xl^i l^if^P = 1- The quantum deple- 
tion is the proportion of particles not in the condensed 
mode (j)"^^'. D 

total number of particles in the system. 



Q 
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particles 



M 



1 



e[o,i], (7) 



and the single-particle von Neumann entropy, 



vN, particles 
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Psp 1 / Psp 



G[0,1], (8) 



1 Nq/N, where TV Tr(pgp*) is the where is the single-particle density matrix: (i|psp|j) = 



(2) 



The pair correlation function, defined as gl 

(f)\V-bjbi), is a measure of the joint probability that a 
particle will be measured at site i and another particle 
will be measured at site j. For visualization, we plot , 
defined as the pair correlation between the center lattice 
site and site k (for the case of an even number of sites 
we take the average of the pair correlation between site 
k and the two centermost sites). 

Out-of-equilibrium time evolution of bosons in an opti- 
cal lattice is an ideal medium in which to study quantum 
entanglement dynamics in a many-body system. We con- 
sider two types of entanglement in characterizing a sys- 
tem of atoms on a Bose-Hubbard lattice: entanglement 
of modes and entanglement of particles. In the former 
case, we can think of each localized Wannier function 
as a quantum system that is, in general, entangled with 
the other M — 1 Wannier functions. Because different 
Wannier functions are spatially distinct, we can refer to 
entanglement between Wannier modes, i.e., sites, as spa- 
tial entanglement or entanglement between sites. On the 
other hand, the particles themselves are interacting with 



iPsp)ij with \i) the lowest-band single-particle Wannier 
state centered at site i. 



IV. QUANTUM SOLITON ENGINEERING 

In this section, we present results obtained by using 
TEBD to simulate the generation and evolution of dark 
solitons in the BHH. To create the dark soliton initial 
conditions, we generalize the methods of density and 
phase engineering originally applied to the continuous 
mean-field NLS [5y] to the discrete mean-field DNLS and 
discrete quantum BHH. Comparisons are made between 
soliton dynamics in the DNLS and BHH. For ah BHH 
simulations in both real and imaginary time, we employ 
a TEBD code optimized to conserve total particle num- 
ber [53. 

Density engineering in a system described by either the 
DNLS or BHH can be accomplished as follows. To mold 
the desired density, one selects an appropriate external 
potential Ck and then finds the ground state of the re- 
sulting model, which we carry out in practice through 
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the form 



FIG. 2: (Color online) Engineered standing soliton profile and 
DNLS dynamics, (a) The density and phase of the engineered 
standing DNLS soliton solution is plotted versus spatial posi- 
tion. Because we select an even number of sites, the two points 
with lowest density straddle the lattice center at x = 0. (b) 
Real time dynamics in the DNLS of the initial condition in 
(a) show a stably propagating density notch with a spray of 
phonons. A plot of the phase (not shown) reveals a tt phase 
drop across the notch for all times. Note that all surface plots 
in this paper have space on the horizonal axis and time on the 
vertical axis. 



imaginary time evolution. On the other hand, phase en- 
gineering is more subtle. In mean- field theory, all atoms 
are assumed to be in the same wave function and hence 
in phase with one another. Thus, imprinting a phase on 
the many-body wave function is conceptually equivalent 
to doing so on a single-particle wave function. Quantum- 
field theory, e.g., the BHH, relaxes the assumption that 
all particles are in-phase. However, we can still imprint a 
phase in the natural way so long as the condensate frac- 
tion is large enough. To do this, we can change the local 
lattice potential to Ck = — ^fe^/^puise for a short time tpuise 
in real time to imprint a phase Ok at site k. In simula- 
tion, it suffices to apply the phase instantaneously. This 
is accomplished by applying the on-site unitary operator 
exp(m^/c) to site /c, where h is the local number operator. 
In the DNLS, applying an instantaneous phase amounts 
to multiplying the DNLS wave function i/j^ by the phase 
factor e^^'' . When applying the phase instantaneously, we 
are assuming that the density is unaltered during phase 
engineering; this requires that tpuise be much less than the 
correlation time h/fi^ where /i is the associated chemical 
potential. The exact forms of the external potential Ck 
used for density engineering and of the imprinted phase 
Ok depend on the problem in question. The forms used 
for soliton engineering are given in the discussion that 
follows. 

We first concentrate on the case of an engineered stand- 
ing soliton in the center of the lattice (see Fig.[T|). A sin- 
gle Gaussian beam is used to engineer the density notch 
and a hyperbolic tangent phase profile is then imprinted 
across the notch. We can dig a single density notch at 
the lattice center, i.e., at x = 0, by performing imaginary 
time relaxation with an external Gaussian potential of 



Ck Vi exp 



(9) 



For the phase imprinting, we use a hyperbolic tangent 
phase profile of the form 



AOi 



tanh 



2xk_ 



1 



(10) 



where in each case Xk/a = k — (M + l)/2 is the position 
of site k with a the lattice constant. 

For a true standing dark soliton, we would have to 
use a Dirac delta function for the potential in Eq. (|9]) 
to effectively create two independent boxes so the wave 
function would heal perfectly to zero at the center of the 
lattice. Thus, for a Gaussian potential of finite amplitude 
and width, a certain fraction of the energy used to form 
the notch goes into the creation of phonons and a certain 
fraction goes into the creation of the soliton [HB]. Also, 
standing solitons have a tt phase drop across the notch 
in the form of a perfect step function; hence, ideally we 
should take AOi = tt and aei 0. Experimentally, the 
parameter a^ii, i.e., the width of the phase profile, rep- 
resents diffraction-limited fall-off of the far-detuned laser 
beam used to raise half the lattice [s^. If agi is on the 
order of the soliton healing length, the engineered soliton 
will drift to the left. We ignore this effect in simulation 
by taking a sufficiently small value of agi- Note, however, 
that stationary solitons can still be created with signifi- 
cant diffraction-limited fall-off by choosing a value of 
slightly greater than tt to create a second shallow soliton 
to carry away the momentum of the first [56l ]. 

Density and phase engineering using the forms given 
by Eqs. ([9|) and (p!Q|) can be applied to either the DNLS 
or the BHH for soliton creation. We consider both cases 
in this section. A schematic of the procedure used to 
quantum engineer a stationary dark soliton is depicted 
in Fig. [TJ 

Let us now consider a single simulation of soliton en- 
gineering in both the mean-field (DNLS) and quantum 
many-body (BHH) pictures. We choose to work on a 
lattice of M = 30 sites with = 30 particles so that 
the filling factor u = N/M = 1. The effective interac- 
tion parameter is taken to be vU / J = 0.35. This is an 
important parameter because, regardless of the values of 
V and U/J separately, the DNLS statics and dynamics 
are equivalent, up to a normalization of the order pa- 
rameter for fixed vU / J at a given M. The relevant 
parameters used to engineering the soliton are Vi/J = 1, 
(yyijob = 1, AO = TT^ (Tei/a = 0.1. The reason that we 
select an even number of sites is that it allows us to set a 
reasonable amplitude Vi for the Gaussian potential. For 
an odd number of sites, Vi would have to be increased 
substantially in order for a node to be created at x = 0, 
a requirement for the dark soliton to be stationary. Not 
increasing Vi for an odd number of sites results in an in- 
stability of the soliton since the wave function does not 
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FIG. 3: (Color online) Average number and number fluctua- 
tions for engineered standing quantum soliton. (a) The aver- 
age particle density during quantum evolution indicates that 
the dark soliton fills in over time; {hk) can be thought of as 
the sum of all natural orbitals weighted by their respective 
occupation numbers, cf. Fig. 2] (b) The normalized number 
variance measuring the relative on-site number fluctuations is 
greatest as the soliton fills in. 



vanish at the origin, an effect we have observed. The soli- 
tons created with this method with an even number of 
sites are similar to the "5 modes" analyzed in Ref. [37l |. 
but without the site-to-site staggered sign structure. 

Figure [21(a) depicts the soliton profile initial condition 
obtained in the DNLS after imaginary time evolution 
with the external potential (|9]) and application of the 
phase imprint (fTQ|) : the subsequent real time dynamics 
are shown in Fig. [21(b). We only show the density in 
Fig. [2Kb), while, as expected, the phase grows approxi- 
mately linearly in time with a tt phase drop across the 
notch for all times. Due to our use of a finite-size Gaus- 
sian for density engineering, phonons are also created. 
We now wish to compare these results to a full quantum 
many-body calculation using TEBD. 

To this end, we take the same parameters as used for 
the DNLS and perform imaginary time propagation in 
TEBD with the same Gaussian potential used to produce 
the density structure in Fig.[2Ka). Retaining x — 120 ba- 
sis sets and allowing up to (i — 1 = 7 particles per site, 
we are able to calculate a ground state well-converged 
density-engineered ground state via TEBD imaginary 
time evolution. The quantum depletion in this ground 
state is only D ^ 11%. We then remove the Gaussian 
potential and model the application of an instantaneous 
phase imprint by applying to each site k the on-site uni- 
tary ex.p{ihOk), where Ok is the same narrow tanh phase 
profile shown in Fig. [l^a). The resulting quantum dy- 
namics are shown in Figs.[3l[H and[5l 

The evolution of the average particle number den- 
sity and normalized number variance can be seen in 
Fig. [3l We see that, in contrast to the DNLS dynam- 
ics of Fig. [21(b), time evolution according to the BHH 
causes the soliton to fill in. This is qualitatively simi- 
lar to the phenomena we have observed using mean-field 
DNLS soliton solutions as initial conditions for quantum 
evolution [13, S^. The normalized variance plot indi- 
cates that the relative number fluctuations are largest in 
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FIG. 4: (Color online) Natural orbital dynamics for engineered 
standing quantum soliton. The densities of the four most 
highly occupied natural orbitals are plotted in (a)-(d) in or- 
der of decreasing occupation Nj. Population of the mode 0^^^* 
is most responsible for the soliton decay shown in Fig. [3ja). 
Here, we show the natural orbital densities l^^"^'*!^ multiplied 
by their respective average occupation numbers Nj . 



the regions of space-time where the soliton is filling in. 

The filling-in phenomena can be attributed to allowed 
two-body scattering processes that deplete particles out 
of the original solitonic condensate wave function into 
modes with nonvanishing density at the lattice center. 
To see this explicitly, we plot the densities of the four 
most highly occupied natural orbitals in Fig. [H We see 
in Fig. [3Ka) that the condensate wave function ^[.^^ re- 
tains its solitonic form throughout time evolution but 
that occupation of higher order modes causes the soliton 
to fill in over time. In the mean-field DNLS theory, the 
symmetry of the initial condition is preserved through- 
out time evolution in a manner very similar to that of the 
single-particle Schrodinger equation. Hence, for an initial 
condition in the form of a dark soliton, the antisymme- 
try of the condensate wave function is preserved during 
mean-field evolution: recall the persistent tt phase drop 
across the notch in the Fig. [21(b). Because mean- field the- 
ory assumes all bosons to be in the same configuration 
for all time, the soliton is a stable configuration of the 
DNLS. 

However, as we have shown, a full quantum treatment 
such as simulation of the BHH predicts depletion out of 
the solitonic condensate wave function into modes that 
effectively fill in the notch. Because derivation of the 
BHH takes into account only two-body scattering pro- 
cesses, it must be these processes that are responsible for 
depletion into non-solitonic orbitals. The symmetry of 
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the many-body wave function is preserved regardless of 
the model chosen. Therefore, the mechanism most re- 
sponsible for soliton graying should involve two atoms 
depleting from the antisymmetric condensate wave func- 
tion into symmetric modes, the lowest energy of which 
is shown in Fig. [H^b) and has significant density in the 
region of the soliton notch. The other main scattering 
process involves two atoms previously depleted into a 
symmetric mode scattering into an antisymmetric mode 
of higher order than the soliton mode. Both of these 
processes retain the symmetry of the many-body wave 
function. Although the population of (j)]^^ is most re- 
sponsible for soliton dissipation, higher order modes also 
contribute, e.g., (j))^ in Fig. [D^c). 

The time dependence of the phase of the natural or- 
bitals is difficult to show graphically because at each time 
step an independent diagonalization of the single-particle 
density matrix is performed, and the diagonalization rou- 
tine arbitrarily assigns a global phase to the eigenvectors. 
However, the symmetry of the natural orbit als can still 
be discerned from these plots regardless of their unsight- 
liness. For completeness, we quote the symmetries of the 
lowest four modes as follows: 0^^^ is an antisymmetric 

mode, (j)"^^ and (j)"^^ are both symmetric, and 0^^^ is an- 
tisymmetric. 

Similar conclusions are reached in Ref. [41] where the 
authors use a three-mode approximation to model the 
antisymmetric to symmetric depletion of a soliton in a 
continuous trap geometry. Our treatment has the ad- 
vantage that bosons are allowed to occupy any of the M 
allowable natural modes of the lattice so that both de- 
pletion processes can occur. Our method also allows us 
to directly probe the quantum many-body nature of the 
system using, for example, the pair correlation function 
and entanglement measures as depicted in Fig. [5l 

By calculating the average local density as in Fig.[3l[a), 
we have only shown that the dark soliton fills in on aver- 
age over many density measurements. Thus, we cannot 
yet say what would happen during a single experiment: 
does the soliton in fact fill in or does it wander randomly 
to the left and right without filling in at all? We can 
address this question by examining the pair correlation 
function g^^^ defined in Sec. Illli If the soliton were wan- 
dering, this measure would vanish at values Xk ^ ^ even 
during the the later times of evolution when {fik) indi- 
cates that the soliton has filled in. This is because is 
a measure of the joint probability for measuring particles 
at the origin (the initial position of the soliton) and at 
site k. But, for a wandering dark soliton, if a particle is 
measured at the origin, there would be another site, 
say, with vanishingly small density due to the presence 
of the displaced soliton. This would imply g)^^ to be 
approximately zero. In Fig. [51(a), we depict the time de- 
pendence of the pair correlation function and show that, 
in fact, there is not some /cq for which g)^^ vanishes even 
during later times. This result is one strong piece of 
evidence that our study is a full quantum entangled dy- 




FIG. 5: (Color online) Quantum measures for engineered 
standing quantum soliton. (a) The pair correlation function 
with respect to the center lattice site tells us that the soliton 
is in fact filling in during a single experiment. The initial 
condition has substantial spatial entanglement as indicated 
by (b) the local von Neumann entropy and (c) the average 
local impurity. At short times, before the soliton has filled in, 
the entropy in the region of the soliton is lower than at other 
points in the lattice, but this notch in entropy fills in over 
time. In (c) and (d), we see that both quantum depletion and 
particle entanglement grow in time. 

namical study of dark solitons; we present yet another 
piece of evidence in Fig. [3 when we discuss a collision 
between two dark solitons. 

In Fig. O we also depict the dynamics of the entangle- 
ment measures defined at the end of Sec. Hill The local 
entropy plot in Fig. [Sfb) shows that the soliton engineer- 
ing causes the local von Neumann entropy to decrease in 
the region of the density notch. Time evolution causes 
this notch in entropy to dissipate. In contrast to our re- 
sults pertaining to quantum evolution of mean-field soli- 
tons discussed in Sec. [Viand Ref. [13], both the local von 
Neumann entropy and average local impurity have finite 
value at t = 0; Qmodes remains approximately constant in 
time although it grows slightly while the soliton is filling 
in. Panels (c) and (d) of Fig. [5l depict this time depen- 
dence of Qmodes as wcll as growth in quantum depletion 
single-particle impurity Qparticies, and single-particle 
von Neumann entropy ^vN, particles- 

The behavior presented above for the quantum evo- 
lution of a density- and phase-engineered standing dark 
soliton is general. The dependence of the dynamics on 
parameter choices can be summarized as follows. Increas- 
ing the effective interaction strength vU / J at fixed filling 
V and lattice size M causes the soliton to decay at a faster 
rate, i.e., the soliton lifetime decreases. Also, increasing 
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FIG. 6: (Color online) Engineered colliding soliton profile and 
DNLS dynamics, (a) As in Fig. [2l we show the DNLS initial 
condition, but this time we plot the state that results in a 
symmetric soliton collision during time evolution, (b) The 
resulting dynamics according to the DNLS reveal an elastic 
collision between the two dark solitons. 



the filling factor at a fixed lattice size, but keeping the 
effective mean-field interaction uU / J fixed, causes an in- 
crease in the soliton lifetime. Since the effective strength 
of quantum fiuctuations is quantified by U/ {2vJ) in the 
ID BHH [sil, [5^, these trends clearly indicate that the 
soliton decay shown above is due to quantum fiuctua- 
tions. A plot of these trends can be found in Fig. 2 of 
Ref. [13] where we use the above method of dark soliton 
engineering in the BHH to compute the soliton lifetime 
versus vU / J at fixed filling factors v = 0.5, 1.0, 1.5,2.0. 

One can apply the very same soliton engineering meth- 
ods used above to create initial conditions in both the 
DNLS and BHH that result in a collision between two 
dark solitons traveling at equal-and-opposite velocities. 
In this case, we use a potential of the form 











jexp 


I 2al2 \ 


+ exp 


[ J 



(11) 



for the density engineering. This potential has the effect 
of creating two identical density notches at positions x = 
±b with depths and widths controlled by V2 and ay 2, 
respectively. We then imprint an instantaneous phase of 
the form 



Ok = A6>2{ - - tanh 



2{xk + b) 



tanh 



0-02 

2{xk - b) 



(^92 



1}- 



(12) 



The speed of the solitons is controlled by the phase drop 
A^2, and phonon generation is minimized by appropri- 
ately tuning the width (Jv2 of the phase profiles to the 
soliton depth as determined by V2 in the density engi- 
neering stage. 

In Fig. [6l we create a soliton-soliton collision in the 
DNLS. The initial condition obtained after density and 
phase engineering is shown in Fig. [61(a), while the sub- 
sequent evolution of the density is shown in Fig. [61(b). 
For parameters, we take vU / J = 0.35 at filling v = 1 
with M = 30 lattice sites and F2/J = 0.4, (Jv2/ci = 1, 
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FIG. 7: (Color online) Quantum-induced inelasticity for quan- 
tum engineered dark solitons. A collision between two engi- 
neered dark quantum solitons is shown for fixed vU / J = 0.35 
at filling factors (a) 1, (b) 0.7, (c) 0.4, and (d) 0.2. The av- 
erage particle number is plotted in each case. The collision 
elasticity decreases with decreased filling due to an increase 
in the growth rate of quantum fluctuations for lower particle 
densities. 



b/a = 6, A6>2 = 0.3 TT, and (792/ cl = 2. The collision is 
elastic, i.e., the speeds of the solitons are unchanged as a 
result of the collision. The forward in time (backward in 
space) trajectory shifts of the solitons due to the collision 
are barely noticeable by eye [60] . We show the results of 
the analogous TEBD simulation in Fig.[7Ka) with numer- 
ical parameters x = 80 and d = 7. For this filling factor 
(i^ = 1), the average particle number density {fik) closely 
follows the DNLS wave function density |V^/cp, especially 
at times before and during the collision. 

However, as occurs for the standing soliton analyzed 
above, there exist two-body scattering processes which 
deplete the system and give rise to entangled quantum 
dynamics. The depletion processes which preserve the 
overall symmetry of the many-body wave function in- 
clude two atoms scattering from a symmetric mode into 
an antisymmetric mode and two atoms scattering from a 
symmetric mode into another symmetric mode of higher 
order. Since antisymmetric modes must vanish at x = 0, 
occupation of such modes will affect a collision of two soli- 
tons that occurs at the center of the lattice. In order to 
see this effect, we can lower the filling factor u while keep- 
ing vU I J fixed in order to increase the rate of quantum 
fiuctuations. Then there is a significant amount of deple- 
tion before or during the collision. In panels (b), (c), and 
(d) of Fig. [3 we show the quantum evolution with the 
same soliton engineering parameters and vU /J = 0.35 
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but with V = 0.7,0.4,0.2, respectively. Because uU/J is 
held fixed, the initial density is the same in each case up 
to a scale factor, but the different growth rates of deple- 
tion give rise to different density profiles at later times. 
In fact, because the mode into which most atoms deplete 
is antisymmetric and vanishes at x = 0, tuning the time 
scale of quantum fluctuations to occur near the time of 
collision effectively induces an inelasticity in the collision. 
This can clearly be seen by comparing the u = 1 case to 
the u = 0.2 case in Fig. [3 We have observed qualitatively 
similar results when considering DNLS mean-field initial 
conditions for subsequent BHH quantum evolution of a 
pair of colliding solitons [13, • 

This result presents another strong piece of evidence 
that a given experiment will not reveal the solitons drift- 
ing to the left or right randomly without changing shape, 
i.e., the solitons are not simply spreading out or diffusing 
without dissipating: if the solitons were spreading out, 
they would not appear to stick together on average over 
many density measurements as in Fig. [3 Hence, the phe- 
nomenon of quantum-induced inelasticity complements 
the result regarding filling-in of the pair correlation func- 
tion in the standing soliton case [see Fig. [51(a)]: when a 
fully quantum many-body simulation of the dark soliton 
is performed, it predicts that the soliton gets depleted 
and fills in, even during a single experiment. 



V. QUANTUM EVOLUTION OF MEAN-FIELD 
SOLITONS 




FIG. 8: (Color online) Natural orbital dynamics during quan- 
tum evolution of a standing mean-field soliton. (a) The con- 
densate wave function and the (b) second, (c) third, and (d) 
fourth most highly occupied natural orbitals during quantum 
evolution of a mean- field dark soliton. Note the different color 
bars in each case. The discontinuities are due to a swapping of 
natural orbital occupation numbers at specific times during 
the evolution as we indicate with dashed-dotted horizontal 
lines. 



In the above analysis, we used density and phase en- 
gineering to generate initial conditions in the BHH that 
resemble dark solitons. However, the dark soliton is the 
first excited stationary state of the DNLS. In this sec- 
tion, we consider what happens when such a mean-field 
soliton state is injected directly into the quantum BHH. 
Two of us have analyzed this case in two other works 
as well [13, Eil, so we will only treat it briefly. Here 
we present new results in terms of Pegg-Barnett phase 
averages, as compared to previous work. 

To summarize the method, we first perform con- 
strained imaginary time relaxation on the DNLS with an 
initial condition in the form of a linear function through 
the center of the lattice. We run the imaginary time evo- 
lution until converged to the fundamental dark soliton 
solution of the DNLS with a notch centered at x = 0. 
We then use a product of on-site coherent states, cf. Eq. 
([2|), truncated at d number states to build a Fock space 
representation of the DNLS soliton solution. The matrix 
product state representation of such a product state is 
trivial, so performing subsequent evolution in TEBD is 
straightforward [46]. It is important to note that in this 
case we must use a non-number-conserving TEBD rou- 
tine because the initial mean- field states themselves, i.e., 
truncated coherent states, do not conserve total particle 
number due to DNLS mean-field theory assuming that 
the U(l) symmetry present in the BHH is spontaneously 



broken. 

In this case, as in Sec lIV) we observe the initial soliton 
configuration to fill in with depleted atoms over time [l^ • 
The analysis of the depletion mechanism is very similar 
to that presented above. In Fig. [HI we show the quan- 
tum dynamics of the four most highly occupied natural 
orbitals, cf. Fig. [H for the analogous plot for density and 
phase engineered initial conditions, using TEBD with pa- 
rameters uU/J = 0.35, = 1, M = 31, X = 50, d = 7. 
The quantity u = A^dnls/^ is now the average filling 
factor with A^dnls being the norm of the initial DNLS 
soliton wave function. We will refer to the plots in Fig. 
[8] when conducting a BDG analysis of the initial DNLS 
dark soliton state in Sec. WH 

Of course, a density notch is not all that characterizes 
a dark soliton; the phase drop across the notch is equally 
important. In fact, in the context of the mean- field GP 
equation, it is exactly this phase drop that stabilizes a 
stationary dark soliton and keeps particles from filling 
it in. Thus, loosely speaking, we can say that quantum 
many-body effects are "stronger" than this phase drop 
in that these effects do cause the soliton to dissipate, as 
we have well shown. In Fig. [9l we consider two different 
measures of the phase of the soliton during quantum evo- 
lution: (1) the average of the Pegg-Barnett phase oper- 
ator, (6>/e), and (2) the argument of the order parameter. 
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FIG. 9: (Color online) Phase measures during quantum evo- 
lution of a standing mean-field soliton. (a) The average of 
the Pegg-Barnett phase operator and (b) the phase of the 
order parameter during quantum evolution of a mean-field 
dark soliton. The phase drop across the notch according to 
the Pegg-Barnett definition decays as the soliton depletes, 
whereas the order parameter phase maintains a tt phase drop 
across x = for all times. 



arg((6/e)). We can also calculate the phase of the nat- 
ural orbitals, i.e., arg((/)^-^^). In fact, the symmetries of 
the natural orbitals were alluded to in the discussion of 
the soliton depletion mechanisms in Sec. lIV| but as men- 
tioned there, it is not insightful to plot the quantities 
arg(0^^^) over time. 

We see in Fig.[9l^a) that at time t = the Pegg-Barnett 
phase operator initially exhibits a phase drop across the 
notch approximately equal to tt. However, when parti- 
cles deplete out of the antisymmetric soliton mode into 
modes of both even and odd symmetry, the phase drop 
across the notch according to the Pegg-Barnett defini- 
tion decreases. For times t > when multiple single- 
particle orbitals are occupied, the phase of the system 
in the single-particle wave function sense is ill-defined 
even though the Pegg-Barnett definition remains valid. 
Of course. Ok being a local observable, this is an effect 
that manifests itself only over many measurements. It 
is also important to state that the on-site Pegg-Barnett 
phase operator as we have defined it requires a breaking 
of the U(l) symmetry, i.e., non-conservation of particle 
number: the only nontrivial part (the second term) of 
Ok in Eq. (|4]) will clearly vanish when averaged over a 
many-body state with conserved particle number. 

In Fig.lHJb), we show the space-time dependence of the 
phase of the on-site order parameter. As the soliton de- 
pletes, the norm of the order parameter, |(^/c)P, decays 
in the tails of the soliton [17]. However, as clearly shown 
in Fig. mb), the tt phase drop of the order parameter 
across the notch persists in time. Next, we seek to an- 
alyze the initial DNLS soliton configurations considered 
in this section within the BDG framework. 



VI. BOGOLIUBOV ANALYSIS 

In this section, we use the Bogoliubov prescription to 
perform a linear stability analysis of the stationary state 
DNLS soliton solutions whose quantum dynamics were 
analyzed in Sec. W\ and Ref. [l^. We make an explicit 
comparison of the spectrum of elementary excitations 
predicted by the Bogoliubov theory versus the depletion 
of the soliton over time as predicted by quasiexact simu- 
lations of the BHH using TEBD. 

Let us assume a stationary state of the DNLS so 
that the amplitude at each site k evolves according to 
^^^-zfit/h^ where /i is the associated chemical potential. 
We then perturb this initial condition and look for solu- 
tions to the DNLS of the form 



• V 



(i) (i) 
where u]^^ and vj:^ 



(13) 

are the amplitudes of the j^^ nor- 
mal mode and coj is the associated mode frequency. In- 
sertion of Eq. (fT3|) into the DNLS results in the tight- 
binding form of the BDG equations which we solve nu- 
merically. The diagonalization routine used to solve the 
BDG equations results in 2M eigenvectors and eigen- 
values. The norm of a given mode j is defined as 
this motivates the definition of 
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E 

the local mode density at site k to be \u\^ 
For modes with nonzero positive norm, we renormalize 



the amplitudes such that 

If for a given mode the frequency Uj is complex, the 
DNLS solution is said to be dynamically unstable. On 
the other hand, if coj is negative and the norm of the 
mode is positive, the mode is anomalous^ in which case 
the DNLS solution is at a local maximum in the en- 
ergy landscape and is said to be energetically unstable. 
A mode with identically zero frequency corresponds to a 
Goldstone mode associated with the breaking of a contin- 
uous symmetry in the stationary mean-field solution. As 
an example, we always observe in our numerical BDG cal- 
culations the Goldstone mode associated with the calcu- 
lated DNLS soliton solution attaining an arbitrary global 
U(l) phase. 

In a harmonic trap in the Thomas- Fermi limit, the 
dark soliton has an anomalous mode with a frequency 
— 1/\/2 times the trapping frequency [iil, EH; this mode 
is associated with translation of the soliton in the trap. 
In a system with full translational symmetry, this mode 
is a Goldstone mode and thus has zero frequency. In 
the continuum, such a zero frequency (and zero norm) 
translational mode also appears in a box geometry given 
that the healing length is much smaller than the box size 
and the soliton is placed sufficiently far from the hard 
walls [l6| . In our case, we should not expect this mode to 
have exactly zero frequency since the lattice itself breaks 
continuous translational symmetry. Under these guide- 
lines, we now proceed to analyze the solutions of the BDG 



,(i)|2 



„(j)|2l 



1. 



11 



0.25 
0.2 

^0.15 

I 

^ 0.1 

_s_ 

0.05 



(a) 





uU 
J 


= 




vU 
J 


= 0.05 




J 


= 0.20 




. vU 

J 


= 0.35 
















-10 



FIG. 10: (Color online) Anomalous Bogoliuhov mode: Density 
profile and comparison to first depleted mode, (a) The den- 
sity of the anomalous (negative frequency and positive norm) 
mode of the BogoUubov excitation spectrum is plotted ver- 
sus space for effective interactions vU / J — 0,0.05,0.20,0.35 
at M = 31 sites, (b) For the latter case of vlJ j J ^ we show 
the time evolution at time steps tJ jft ^ 25, 30, 35, 40 of the 
depleted natural orbital most responsible for the decay of the 
dark soliton during the initial stages of evolution. We see 
that the Bogoliubov mode in (a) for uU/J = 0.35 does not 
resemble the evolving symmetric mode shown in (b). 



equations for a stationary DNLS soliton background. 

For the case of M = 31 sites as in Sec. IVland Ref. [13], 
for interaction strengths of interest, vU/J < 0.50, we ob- 
serve in our calculations an anomalous mode with small 
negative frequency —0.05 J/ h < cji < 0. As expected, 
this anomalous mode is characterized by a density max- 
ima concentrated in the region of the soliton notch as 
shown in Fig. [TOlf a) for various values of uU / J. We wish 
to compare the density of this anomalous Bogoliubov 
mode to the most highly occupied symmetric natural or- 
bital during full quantum evolution. Figure [TOlfb) depicts 
the density of the symmetric natural orbital with highest 
occupation at times tJ/h ~ 25, 30, 35, 40 during quantum 
evolution according to the BHH of the stationary DNLS 
soliton [see Fig. [8l(a)-(b)]. Strictly speaking, this mode 

is for tJ/h < 33 and for tJ/^ ^ 33 due to the 
ordering of the natural orbitals being defined to coincide 
with the ordering of occupation numbers. Recall that for 
this simulation vU / J = 0.35. It is clear that the sym- 
metric natural orbital in Fig. [TOlf b). the natural orbital 
most responsible for depletion out of the soliton orbital, 
does not closely resemble the negative frequency mode 
of the Bogoliubov theory: this natural orbital, unlike the 
anomalous Bogoliubov mode has significant density away 
from the lattice center. Also, as can be seen in Fig. [HI 
there is also significant contribution from natural orbitals 
of higher order, e.g., , , etc., to the filling-in ef- 
fect. A truncation of the Bogoliubov Hamiltonian at the 
anomalous mode would thus hardly be justified in this 
regime of strong quantum fluctuations. 

The dependence of the frequency of the anomalous Bo- 
goliubov mode versus effective interaction strength vU / J 
for the case of M = 31 sites is depicted as the solid blue 
curve connecting circle points in Fig. [Til We can un- 
derstand the form of this dependence as follows. The 
anomalous mode is associated with the translation of the 
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FIG. 11: (Color online) Frequency of Bogoliuhov mode versus 
interaction strength: Appearance of energetic and dynamical 
instability. For M = 31 and M = 30 sites, we show the real 
and imaginary parts of the calculated negative or complex 
frequencies appearing in the Bogoliubov excitation spectrum. 
For an odd number of sites, e.g., M = 31, the frequency is 
always negative and real, while for an even number of sites, 
e.g., M — 30, the frequency becomes purely imaginary past a 
critical value of vU/J. This behavior is general for any odd 
or even number of sites and is not specific to the particular 
choice of lattice sizes made above. 



soliton in our lattice plus box geometry and in general 
could have nonzero negative frequency due to the lattice 
and box breaking full translational symmetry. For the 
case of vU / J = 0, our "soliton" state is the first excited 
state of the linear tight-binding Schrodinger equation. 
The size of the density notch is thus on the order of the 
system size and the box works to break the translational 
symmetry associated with the density notch position in 
a manner similar to how a harmonic trap breaks this 
symmetry: for vU / J = 0, the translational mode is al- 
ways anomalous. As we increase vU / J, we are decreasing 
the healing length, i.e., the soliton size, and thus restor- 
ing translational symmetry of the soliton inside the box; 
thus, the frequency of the anomalous mode becomes less 
negative. In the continuum, the mode would eventually 
become a zero mode, cf. Ref. [l6j. However, as we con- 
tinue to increase uU / J and decrease the soliton size, the 
breaking of translational symmetry by the lattice begins 
to play a role. Past some critical value of vU/ J at which 
uji is nearly zero, the frequency goes more negative upon 
further increasing uU / J. This picture is general for the 
case of an odd number of sites in which the soliton notch 
is located at a particular site. Recall that for simplicity, 
we are here only considering solitons placed at the lat- 
tice center. Also, as expected, the critical value of uU / J, 
past which coi begins to decrease, decreases as the sys- 
tem size M is increased. This is because the soliton size 
goes like J/uU^ and the larger the lattice size is, the 
more quickly translational symmetry within the box is 
restored as vU/J is increased from zero. 

On the other hand, for the case of an even number of 
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sites in which the notch of the sohtonic stationary state 
of the DNLS is located between sites, the behavior is 
quahtatively different from that described above for an 
odd number of sites. As shown in Fig. [11] for M = 30, 
the negative frequency mode goes to zero at some crit- 
ical value of vU / J. Now, rather than becoming more 
negative past this critical value as for M = 31, the mode 
instead becomes complex, thus indicating the onset of 
a discreteness-induced dynamical instability. The emer- 
gence of such a dynamical instability has been pointed 
out also for a dark soliton state in a double-well poten- 
tial [H, [13 , which corresponds to the limit of the small 
even number of sites, M = 2. However, we have checked 
that quantum evolution of DNLS dark solitons for an 
even number of sites is qualitatively identical to the be- 
havior depicted in Sec. [VI and Ref. [13] in which an odd 
number of sites was assumed. We could also use the imag- 
inary part of the complex frequency mode to estimate the 
timescale over which dynamical instability would mani- 
fest itself for the soliton engineering calculations in Sec. 
IIV| which for practical reasons were performed with an 
even number of sites. Using the data in Fig. [Til we 
estimate that the timescale of dynamical instability for 
vU / J = 0.35 and M = 30 sites is three to four times 
longer than the longest times shown in the characteristic 
simulation of Sec. [TVl 

We have also considered higher order Bogoliubov 
modes, e.g., {u^^\v^^'^)^ {u^^\v^^'^)^ etc., for both the soli- 
ton and DNLS ground state background configurations. 
According to our calculations, the lowest natural orbitals 
obtained using TEBD do not match the low energy Bo- 
goliubov modes of the formed dark soliton state, although 

(2) 

the second and third depleted natural orbitals, i.e., 0^ 

and 0^ , compare more favorably to the second and third 
lowest energy Bogoliubov modes calculated with respect 
to the DNLS ground state Stih, as in Fig. [lO^a)- 

(b), direct comparison is made difficult because the nat- 
ural orbitals themselves are time-dependent objects. Ex- 
actly how the noncondensed atoms distribute themselves 
among the single-particle natural orbitals cannot be ac- 
curately predicted by calculating the Bogoliubov quasi- 
particle spectrum of DNLS stationary states. Finally, 
such a static Bogoliubov analysis as we have presented 
in this section cannot strictly be applied to the engi- 
neered dark solitons of Sec. IIVI because the density and 
phase engineered initial conditions encountered in that 
case are not stationary states of even the DNLS. It would 
in principle be possible to carry out time-dependent BDG 
calculations in the spirit of Ref. [43], treating the soli- 
ton engineering as part of the dynamical process. How- 
ever, the substantial amount of quantum depletion and 
quantum entanglement present as predicted by our full 
quantum TEBD simulations suggests that such a theory 
based on noninteracting Bogoliubov quasiparticles would 
inevitably be inapplicable. On a related note, breakdown 
of BDG theory when applied to ground state properties of 
the ID BHH is known to occur for interaction strengths 



as low as U/J> 0.2 



VII. CONCLUSIONS 

We have shown that dark solitons decay in the presence 
of strong quantum fluctuations induced by an optical lat- 
tice, even in the superfluid regime of the Bose- Hubbard 
Hamiltonian. Our work is complementary to past mainly 
perturbative studies of quantum fluctuations in which it 
was found that dark solitons delocalize but do not decay. 

Specifically, we used time-evolving block decimation 
to first engineer one or more dark solitons in an experi- 
mentally accessible manner through manipulation of the 
lattice potential itself. We then followed the entangled 
dynamics of this decay by evaluating a large suite of 
quantum measures, including several entanglement mea- 
sures, moments of number and Pegg-Barnett phase oper- 
ators, quantum depletion, and second-order correlations 
in the form of We showed that not only the mean 
field but also the Boguliubov-de-Gennes equations break 
down, even in the superfluid regime of the phase dia- 
gram. We made an explicit comparison of BDG analy- 
sis and our entangled dynamical simulations, and found 
that "depletion modes," i.e., higher natural orbitals in 
our more complete theory, do not in general resemble 
the BDG modes. 

A main aim of this paper is to make previous work 
presented in Ref. [13] relevant to realistic near-term ex- 
periments, as well as to compare our work in a more 
explicit way to past BDG-based studies [13, [l5|, \Tdl\ . 
The density and phase-engineering simulations of Sec. IIVI 
can be realized in the lab with experimental parame- 
ters as follows. The Bose-Hubbard effective interaction 
parameter vU / J = 0.35 at filling v = 1 corresponds 
to ^'^Rb atoms with an s-wave scattering length tuned 
via Feshbach resonance to as = 1.0 nm in an optical 
lattice with longitudinal and transverse lattice heights 
Vq = Eji and Vo± = 25 created with lasers of wave- 
length A = 850 nm, where En is the recoil energy of 
the lattice. The Gaussian beam used for density engi- 
neering then has height Vi = J ^ 0.17 En and width 
avi = a = 425 nm. The width of the fall-off of the 
tanh phase profile is aei = 0.1a = 42.5 nm. Increasing 
the width of the phase profile to the order of the soliton 
healing length will cause the soliton to drift to the left, 
as we have observed in simulations. 

In future studies we would like to consider the contin- 
uum limit as suggested by Fleischhauer [13, [65[. Then 
there is an interesting interplay between holding uU / J 
constant in order to keep a consistent ratio of nonlinear- 
ity to dispersion in the mean-field theory, while at the 
same time holding U/vJ constant in order to reach the 
continuum. This would allow us to make explicit ex- 
perimental predictions for the Heidelberg and Hamburg 
groups; at present our simulations remain in the weakly 
discretized limit and cannot be clearly and rigorously ex- 
trapolated to recent experiments. On the other hand, our 
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predictions are relevant for the specific context we have 
considered, namely, a ID lattice in the strongly quantum 
limit of unity filling or less. Such systems have already 
been achieved in ground-breaking experiments [6^, [67l |: 
in this context solitons would be engineered in an array 
of ID tubes. Then one runs O(IOO^) dark soliton experi- 
ments simultaneously. Such high statistics can be useful 
in extracting ^^^^ from noise correlations [68]. 

We acknowledge use of the time-evolving block decima- 
tion algorithm from the TEBD open source project [6^ 
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